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Abstract. In this paper, we investigate the numerous parameters choi- 
ces for linear lattice Boltzmann schemes according to the definition of the 
isotropic order given in [3]. This property — written in a general frame- 
work including all of the D^Q^ schemes — can be read through a group 
operation. It implies some relations on the parameters of the scheme 
(equilibrium states and relaxation times) that give rigorous methodology 
to select them according to the desired order of isotropy. For acoustic ap- 
plications in two spaces dimensions (namely D2Q9 and D2Q13 schemes) 
this methodology is used to propose a full description of the sets of pa- 
rameters that involve isotropy of order m (m G {1, 2, 3, 5} for D2Q9 and 
m G {1, 2} for D2Q13). We then propose numerical illustrations for the 
D2Q9 scheme. 



Introduction 

Lattice Boltzmann schemes [3 W\\ [8] are numerically very interesting be- 
cause of their efficiency and then represent a promising field in computation 
fluid dynamics [SllI^- However, since the directions of the lattice are privi- 
leged, these schemes are a priori not isotropic even though they are used to 
solve isotropic phenomena. All the more, [TD] mentions that isotropy also 
improves the lattice Boltzmann scheme in the sense of the stability. Thus, 
the isotropy properties of the ^dQq schemes (where d represents the space 
dimension and q the number of discrete velocities) are the aim of numerous 
previous works though isotropy is not defined in the same way: for instance, 
isotropy of the stress tensor required by fluid equations in [5], isotropy by 
solving an eigenvalues problem in 

As it is explained in [3], because of the criss-cross pattern of the scheme, it 
is natural to look for parameters that give the same behavior on each axes. 
By extend, isotropy can be related to the invariance of the spatial frame 
by all the special orthogonal transformations. Space transformations are 
actually considered in [131 El SI for discrete cases: isometrics of the Bravais 
lattice in [13] , discrete rotations in [TH U] . An abstract way is considered 
in [3] in order to obtain a systematic procedure to construct higher-order 
BGK models. 

In [3], we give a rigorous definition for a D^Q^ scheme to be isotropic. 
This definition is based on an invariance of the equivalent equations with 
regard to all the special orthogonal transformations of the space and can 
be used for every scheme. Moreover, it gives a procedure to precise the 
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parameters (namely equilibrium states and relaxation times) of linear and 
orthogonalized schemes in order to improve them in the sense of the isotropy. 
This methodology finds some set of coefficients again that are already known 
thanks to the experiment. 

In the first Section, we recall some notations about the linear D^Qq scheme 
defined with orthogonalized moments and then the definition of the equiva- 
lent equations that come from a Taylor expansion of the lattice Boltzmann 
scheme [8j. 

In the second Section, we briefly recall the definition of the isotropy given 
in [3]. A DrfQq scheme is said isotropic at order m if the equivalent equa- 
tions of order m are isotropic. This definition is written thanks to a group 
operation and gives an easy systematic procedure to precise some of the 
parameters of the linear lattice Boltzmann scheme. The end of this Section 
consists in the calculation of the degrees of freedom that we have to take 
into account for every scheme. 

In the third Section, we focus on acoustic applications using the D2Q9 
and D2Q13 schemes. We first recall the definition of this particular schemes 
(moments and equilibrium states) and the methodology of the formal com- 
putation. Then we give a total description of the sets of parameters that 
involve isotropy of order m {m from 1 to 5 for the D2Q9 scheme and to 2 
for the D2Q13 scheme). Concerning the D2Q13 scheme, we also propose a 
selected choice of sets of parameters that involve isotropy at third order. 

In the last Section, preliminary numerical results are discussed to illus- 
trate the different orders of isotropy for the D2Q9 scheme. 



1. Lattice Boltzmann method and equivalent equations 

In this Section, we first recall some notations about the lattice Boltzmann 
method and more precisely we focus on a linear LB scheme with orthogo- 
nalized moments. We then give the formal equivalent equations of order m 
that is the partial differential equations with which the scheme is consistent 
at order m, for all integer m. 



1.1. Notations for the lattice Boltzmann method 

We use the notation proposed by d'Humieres in [7]: we consider a regular 
lattice in d dimensions C with typical mesh size Ax. The time step At is 
determined thanks to the velocity scale A by the relation: 

Ax 

At = - 

For the g-velocities scheme denoted by D^Qg, we introduce V = (^^j)o<j<g-i 
the set of q velocities and we assume that for each node x of £, and each Vj 
in V , the point x + Vj/S.t is also a node of the lattice C. 

The aim of the D^^Q^ scheme is to precise the particle distribution / = 
(/j')o<j<(j-i foi' 3; G i2 and discrete values of time t by solving the following 
PDEs: " 

(1) ^+^rV/, = -;^(/,-/f), 0<j<g-l, 
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where f^'^ describes the distribution fj at the equihbrium and tj is the 
relaxation time (related on fj). 

The interest of the D^Q^ scheme lives in the method of resolution of ([1]): 
we proceed in two steps like in a splitting method as explained below. 

• First, we are interested in the relaxation /—)•/* step that consists 
in solving 

df* 1 

(2) ^ = --(/^.-/;q), 0<j<g-l. 

In order to solve we first consider the moments at equilibrium. 
These moments are divided into two types: the ones that are con- 
served at equilibrium denoted by W G and the ones that are 
not conserved at equilibrium denoted by 1^ S M.'^~^ . More precisely, 
because of the acoustic applications, the moments W and Y are lin- 
ear combinations of the distribution functions / (i.e. there exists an 
invertible matrix M = (-Mjj)o<jj<g-i such that {W^,Y^f = M f). 
The size and the definition of the matrix M depend on the scheme, 
however for every D^Q^ scheme the N first lines of M are the same. 
In fact we have W = {p,q^)^ where the total density p and the 
momentum q are given by 

q-l q-l 

Wo = P = ^fj and VF^ = q„ = ^ Vjfj, 1 < a < d. 

j=0 j=0 

Now, we define the moments after the relaxation thanks to the rela- 
tions: 

W,^ = Wk, 0<k<N-l 
l^fc* = Yk +SkiY^'^-Yk), N<k<q-1, 

where, for N < k < q—1, Sk is related to all the tj, N < j < q—1 
and Y^^ is the moment at equilibrium. For stability reasons, we 
choose Sk such that < < 2, iV < A; < (7— 1. In order to simplify 
the calculations, we introduce the reals a^, N < k < q—1, by 

1 1 

(4) CTfc = -, CTfc > 0. 

As the scheme is supposed to be linear, the moments at equilibrium 
y^'^ linearly depend on the conserved moments W through a matrix 
denoted by i?, so that equation ([3]) reads: 

^*)=j(^), with J=[^^l Id,_°-S 

where S is the diagonal matrix of the relaxation times Sk, N < k < 
q—1. Finally, we reconstruct the distribution after the relaxation 
with/* = M-^ {W*^,Y*^f. 



(3) 



Second, we solve the transport step 
dt 



df- 

(5) ^ + • V/, = 0, 0<j<q-l, 
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thanks to the characteristics method: 

(6) fj{x + VjAt,t + At) = f*{x,t), Vx G £, < j < 1. 
Finally a time step of the linear lattice Boltzmann scheme D^Q^ reduces 

to: 

(7) f{x + VjAt,t + At) = M-^JMf{x,t), Vx G £, < j < g- 1. 

1.2. Equivalent equations 

By convention, we denote •* the ith line and ' j the jth column of a tensor 
• . By convention, a Latin letter is an index related to the moments (of size 
A'^) while Greek letter represents an index related to the space dimension (of 
size d). We then introduce the system of A'^ equivalent equations [8] of order 
m with which the linear lattice Boltzmann scheme is consistent at order m 
(that is to say that the rest is ©(At"*)): 

m 

(8) dtW' + ^ V"W = 0, 1 < i < iV, 

n=l 

where S Sn,d,Nj 1 < n < ni, 1 < i < N are tensors of order n + 2 that 
take into account the coefficients of E and 5. The space Sn,d,N is defined as 
the quotient space Sn,d,N = '^"Z ~> where two tensors are equivalent for 
the equivalence relation ~ if the associate partial differential operators are 
the same (using Schwarz's theorem). By convention, the maximal contrac- 
tion operator is defined by, 1 < i < N, 1 < n < m: 

l<j<N 

1<Q!1,--- ,an<d 

These equivalent equations come from a formal calculus explained in [TO] 
and an algorithm that is easy to use is described in [2]. 

2. ISOTROPY CONDITION 

In this Section, we first briefly recall the construction of the definition of 
isotropy for a D^^Q^ scheme given in [3]. The D^^Q^ scheme is characterized 
through the set of its equivalent equations. And as these equations are 
PDEs, it is natural to derive this definition from the definition of an isotropic 
system of PDEs. Consequently the isotropy reads as a group operation and 
involves a set of equations on the coefficients of the matrices E and S. The 
second step of this Section consists in counting the number of degrees of 
freedom given by the isotropy. 

2.1. Isotropy and Algebra 

We first recall the definition of an isotropic system of PDEs: 

Definition 1. The system of PDEs 

M 

(9) dtW' + ^ = 0, 1 < i < iV, 

n=l 
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is said isotropic if it is invariant by special orthogonal transformation of the 
frame. 

For r an special orthogonal transformation of the frame in d spaces di- 
mensions, we define the orthogonal matrix R{r) by 

such that R{r)^^W is the vector of the conserved moments in the new frame. 

Remark 2. In this study, we focus on the group of special orthogonal trans- 
formations in two dimensions that is the group of rotations SO2 (M) • It can 
he parametrized with one real number corresponding to the rotation angle, 
so that r reads 

( cos(0) -sin(e) \ 
V sin(0) cos(0) j ' 

As it is proven in [3], the isotropy property can be read through a group 
operation defined in the following Definition. Furthermore, this definition 
gives a set of equations that are really easy to use and that give some 
relations between the parameters of the scheme. 

Definition 3. Let n be in N*. Then the group operation • *5'0f^(M) x 
Sn,d,N — ^ Sn,d,N is defined for l<i,j<N,l< ai, • • • , a„ < d, 1 < n < m, 
by the relation: 

(10) («>„(r)(„A))f 1'-'"- := 

l</3l,--- ,l3n<d l<k,l<N 

We then obtain a characterization for a PDF to be isotropic (for the proof 
see [3]): 

Proposition 4. Let Q be a PDE of order m. It is isotropic if the tensor n-A 
is a fixed point of 1 < n < m, that is ^>„(r)(„A) = Vr G S'Orf(M), 
1 < n < m. 

Finally, we recall the definition for a lattice Boltzmann scheme to be 
isotropic 0: 

Definition 5. A lattice Boltzmann scheme is said isotropic at order m if 
the system of equivalent equations ^ at order m is isotropic. 
Furthermore, we denote by L]\f{r) := J2i<n<Ni^'n.{'''){nA) — „yl)At'" the 
lack of isotropy at order N for the special orthogonal transformation r. 

2.2. Consequences on the degrees of freedom 

In this Subsection, we count the degrees of freedom that have to be taken 
into account for studying isotropy. The total number of the parameters that 
can be freely chosen is indeed reduced by considering the relations imposed 
by the physics. 

Let us make the count of the equations and of the unknowns: 

• There are iV^ (d+n — 1)1/ {n\{d— 1)1) different components of each 
one corresponding to a relation for isotropy thanks to Proposition |3] 
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• Then, we take into account the equihbrium states, namely the co- 
efficients of the matrix E defined in Section ll.li The matrix E 
is usuaUy chosen sparse due to the experiments and to the physi- 
cal "good sense". It is remarkable, that the isotropy condition for 
the DdQq scheme could give some justifications of certain traditional 
choices of coefficients. However, in this paper, E is chosen full in 
order to give as degrees of freedom as possible, given then N(q — N) 
additional parameters. 

• It remains to consider the relaxations times. Actually, as it is re- 
called for example in [3] , there exists essentially three types of lattice 
Boltzmann schemes according to the values of the relaxation times: 
the Bhatnagar - Gross - Krook (BGK) scheme, the Two Relaxation 
Times (TRT) scheme and the Multiple Relaxation Times (MRT) 
scheme. Since our purpose is precise the relaxation times in order to 
improve the behavior of the scheme in the sense of the isotropy, we 
consider the third case that is the most general (even if it is not the 
most often used). Thus, we consider q — N independent unknown 
parameters on S. 

Finally, as we consider the most general D^Qq scheme in the linear case, 
there are {q — N){N + 1) free parameters and N^{d + n — l)l/{nl{d — 1)!) 
equations given by the isotropy conditions. However, all of these parameters 
do not have to be precise by the isotropy conditions: some are linked with 
some physical properties through equivalent equations [8]. 

• The linear equilibrium state of the energy described by E§pX'^ + 
-Ef^gxA^ + Es^QyX"^ + Ej'qzX^ is such that the density contribution 
Ee is linked with the sound velocity cq. For example, as it is recalled 
in [ni [13 El [TO], we have E^ = 6cg - 4 for the D2Q9 scheme, E^ = 
2Qcl - 28 for the D2Q13 scheme, E^ = STcg/A^ - 30 for the D3Q19 
scheme, and E^ = 3cg — 2 for the D3Q27 scheme. In fact, all of these 
relations are given by the physicals properties that imply the values 

= = Cg, namely the contribution of dp/dx (respectively 
dp/dy) in the conservation equation that describes the evolution of 
the momentum q^ (respectively qy). 

• Furthermore, since there exists a D2Q9 scheme consistent with the 
Navier-Stokes equations (see [8]), a link shall be done between both 
the viscosities of the considered fluid and two of the relaxation times 
(see for example [101 [TT] ) . Namely, if the shear viscosity is denoted 
by ^ and the bulk viscosity by /i, in the case of the D2Q9 scheme, we 
get 

(11) /i = iA Axfj^^^ and C = A AxfTe(| - Cq), 

where Cg is the sound velocity, a^p^^ (respectively as) is related to 
the relaxation time depending on the moment m^^^ (respectively 
m^) thanks to (Q, where each moment is defined below. 

• For acoustic applications, our objective is to conserve the freedom 
of both viscosities even if they are sometimes taken equal. In other 
words, we are looking for a set of coefficients that does not link Sg 
and s^^^. 
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• Moreover, since the moments and m^^^ are of the same order 
(namely 2), we have to take into account a lattice Boltzmann scheme 
with multiple relaxation times for acoustic applications (a Two Re- 
laxations Times scheme leads to a link between and C because 
and s^^^ are equal). 

Finally it remains {q — N)(N + 1) — 3 parameters to be precise. 

3. ISOTROPY FOR ACOUSTIC APPLICATIONS IN TWO DIMENSIONS 

The aim of this section is to investigate the property of isotropy at differ- 
ent orders for two schemes often used for acoustic applications (D2Q9 and 
D2Q13) s-iid in particular to determine all of the possible choices of coeffi- 
cients that yield an isotropic scheme of order m. We first recall the definition 
of the D2Q9 and the D2Q13 schemes through their moments. After specifying 
the degrees of freedom, we explain the used method to determine all of the 
solutions of these huge non linear systems. In the second part of this section, 
we explicate the results order after order for the D2Q9 scheme and in the 
third part, we precise all of the obtained results for the D2Q13 scheme, the 
intricacy of the equations increasing so much with the order that we have 
to restrict ourselves to the second order. 

3.1. Generalities and methodology 

In this subsection, we specify the notations for the D2Q9 and the D2Q13 
schemes and we propose a methodology to solve the very large systems of 
non linear equations that appear when the isotropy is investigate at each 
order. These systems are written in terms of the coefficients of the matrices 
E (describing the equilibrium states) and S (describing the relaxation times) 
and have to be satisfy for every rotation angle 6. 

For both considered lattice Boltzmann schemes (D2Q9 and D2Q13), three 
moments are conserved during the collision step: the density p and the two 
coordinates of the macroscopic momentum and qy. They are defined by 



Next, we introduce six moments that are not conserved during the collision 
step: the kinetic energy m^, the square of the kinetic energy rngj, the coor- 
dinates of the heat fiux m^^ and my:,^ , and two moments of order two m^^^ 
and m^^ . They read 



q-l 

3=0 



q-l q~l 



q-l 



q-l 



j=0 

q-l 



j=0 



q-l 



nit 



771, 



j=0 



j=0 
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In this paper, we consider the hnear D2Q9 scheme obtained after or- 
thogonahzation of these nine moments, with a full matrix E describing the 
equilibrium states and a full diagonal matrix S describing the relaxation 
times that are written 
(12) 



E :-- 



( E§\^ 


El^X 


El^XX 




EtlX^ 


e!ix^ 


E'^xX' 


EtX^ 


EllX^ 




E%\' 


EfyX' 


EfpxxX 


TpQx \ 


EifxxX 




^iPxy'^ 


E^lyXj 



S := Diag(s£ 



■^£2 ' ^iPx 1 ^iPy ) ^'fixx 1 ^'fixy ) 



In order to define the D2Q13 scheme, we have to introduce four additional 
moments: the cube of the kinetic energy meg, two moments of order five 
rrixe^ and mys2, and finally a moment of order four rrixxe- 



m, 



£3 



E 

0<i<g-l 



1|„, |2\3 



2^ 



) /.> 



m 



m 



— ^ ( 2 1 I ) fj ' 

0<i<g-l 



xe2 — ^ '^^ /i' 
o<j<<?-i 



rric 



0<i<g-l 



In this paper, we consider the linear D2Q13 scheme obtained after orthog- 
onalization of these thirteen moments, with a full matrix E describing the 
equilibrium states and a full diagonal matrix S describing the relaxation 
times that are written 



E 



/ E'eX^ 


Ei^X 


EfX \ 


EfpxxX'^ 


E^lxX 


E^xxX 


EfpxyX 


^V>xy'^ 


E'^xyX 


E'^xx:' 


EtX^ 


EtX^ 


^'Py^ 






E^EiX 


Ell^X"" 


rply \ 4 


Ey£2 X^ 


Eye2 A^ 


rply \4 

njy£2A 


E§,X^ 


El^X^ 


EtlX^ 


^?^\ 


EllX^ 


EllX^ 


\ExxeX 


pqx 1,3 


ExxeX'^ / 



(13) 



S := Diag(se,s, 



Vxx 1 ^<fixy ) ^' 



^x ; ' 



'X£2 ' °y'^2 ' ■'£2 ' ''£3 ' ^xxe 



)■ 



We then apply Definition [S] on the D2Q9 scheme in order to describe all 
of the sets of parameters that improve these schemes (in the sense of the 
isotropy). In two space dimensions, the special orthogonal transformations 
r are the well-known rotation matrices, = 3 and m < 5. We first establish 
the degrees of freedom for each schemes at each order. Section [2.21 gives the 
number of equations (this number does not depend on the scheme): 9(n + l) 
where n is the considered order. The explicit calculation of these numbers 
is given in Table [TJ Furthermore, there are 21 unknowns parameters for the 
D2Q9 scheme and 37 for the D2Q13 scheme. 
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order 


1 


2 


3 


4 


polynomial equations 


18 


27 


36 


45 



Table 1. Number of polynomials eqs. for two dimensions 



We remark that all of these equations are polynomial in cos 6 and sin 6, 
where 6 is the angle that parametrizes the rotation. The study of their 
coefficients then gives other equations that only depend on the coefficients 
of S and E. For instance, for the D2Q9 scheme, at first order we get 18 
polynomial equations thanks to Proposition [5] and the formula: 

= E E {R{r)yiiA'fr'^^{R{rr')l 1 < ^, j < 3, 1 < a < 2. 

I<l3<d l<k,l<N 

All of these equations shall be filed as follow: 

• 6 of them are null. 

• 8 of them are of type a cos^ 6 + b cos^ 6 sin 9 + c cos 6 + d sin ^ + e = 0. 
Since the functions 9 1— t- cos^ 9,6^ cos^ 9 sin 9,9^ cos 9,9^ sin 0, 
and 9^1 are free, each of these equations implies 5 additional 
equations: a = h = c = d = e = Q. 

• 4 of them are of type a cos sin + 6sin^ = 0. Since the functions 
9 I— 7- cos 9 sin 9 and 9 1— ?• sin^ 9 are free, each of these equations implies 
2 additional equations: a = 6 = 0. 

The isotropy for the D2Q9 scheme at first order is characterized by 48 
equations: 8 coefficients of the matrix E are then fixed and it remains 
21—8 = 13 parameters for the isotropy at second order. 

Following this methodology order after order, we are able to list in Table[2] 
and [3] for each order the number of equations that we have to take into 
account (namely, equations that do not depend on the angle 9) and the 
number of parameters that we have to precise. 



order 


1 


2 


3 


4 


equations 


48 


78 


116 


169 


parameters 


21 


13 


8 


6 



order 


1 


2 


3 


equations 


48 


78 


148 


parameters 


37 


29 


23 



Table 2. Number 
of eqs. for D2Q9 



Table 3. Number 
of eqs. for D2Q13 



Remark 6. Both numbers of equations for the D2Q9 and D2Q13 schemes 
are the same (the order n being fixed) because the tensors „^ have the same 
size. However, the coefficients of these polynomial equations (in sva.{9) and 
cos{9) ) depend on the geometry of the scheme and as a result the number of 
equations we have to solve is not the same for the two schemes. 

The methodology to investigate the isotropy at order n is then resume in 
three steps. 

(1) compute the equivalent equations of order n (they could be obtained 
using a formal code [10^ [2]) of the lattice Boltzmann scheme de- 
scribed by the matrices E and S taking into account the relations 
obtained at previous orders. 
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(2) write the whole of the equations given by Proposition [U the number 
of equations being given in Table [2] and Table [3l 

(3) precise the relations on the coefficients that have to be satisfied in 
order to have the desired order of isotropy. 

Finally, we introduce a notation that yields to identify the equations: 

Notation 7 (eqal. . .anj). A natural way to identify the equations consists 
in specify the number of the conservation equation and the corresponding 
coefficient of the tensor i^fA, 1 < N < m. Then, we denote by eqial . . .anj 
the equation given thanks to n-^f^^' W^e then denote the coefficient of 

cos'' (0) sin' (0) in the equation eqal. . .anj by eqia. . . anj cosksinl. 

3.2. Results for the D2Q9 scheme 

In this Subsection, the main result on the D2Q9 scheme is given in Propo- 
sition [HI The proof is then subdivided into five Lemmas by proceeding order 
after order. 

Proposition 8. Let L^(r) be the lack of isotropy for the D2Q9 scheme at 
fifth order for the rotation r, then we get the following propositions. 

• The scheme is isotropic at first order, that is L^{r) = C'(At^),Vr E 
502(M), iffEr, El\ E^^^, Ell^, Ell^, E'^^^, E^l^, and E^^ are 
zero. 

• The scheme is isotropic at second order, that is L^{r) = 0{/S.t^),\/ r S 
502 (I^) , iff a is isotropic at first order, E^^ , E^f^ , E^y , E'^y are zero. 



anc 



d Et = Ell = - 4^^.J/(2^^., + ^. 



• The scheme is isotropic at third order, that is L^{r) = 0(At ),Vr G 
502(M), iff the three following properties are satisfied: 

— it is isotropic at second order, 

- <^'P.. = (^vxy (^^d Ell = -1, 

— at least one of these three properties is satisfied 

* 2E§, +4 + 3E§ = 0, Etl = Ell = 0, 

* cT^. = <y^y = l/(12c7^.J, Ell = EZ = 0, 

• The scheme is isotropic at fourth order, that is L^{r) = ©(At^), Vr G 
502(1^); iff the three following properties are satisfied: 

— it is isotropic at third order, 

- 2EP, + 4 + 3E§ = 0,ae = a^^^, El^ = Ell = 0, = = 

- 2 + 3EP = Q orae,=a^^^. 

• The scheme is never isotropic at fifth order, that is to say the prop- 
erty L^{r) = C'(At^),Vr G S02{^), is never true. 

Remark 9. Some of these properties that insure the isotropy of the D2Q9 
scheme are well-known and often used. In particular, the parameters that 
are taken null in order to obtain the isotropy at first and second orders are 
also null for isotropic reasons considering the kinetic solution at equilibrium 
in the continuous environment. However, the results on the third and fourth 
orders are more surprising even if some usual D2Q9 schemes can be seen as 
peculiar cases of those we propose. 
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The proof of Proposition [8] is detailed below thanks to five Lemmas: one 
for each considered order. The followed methodology is given in the previous 
subsection. We have to remark that we continually use the fact that the 



coefficients 



'^£2' '^<Px^ ^V>y^ '^<PxxJ 



, and aip^y are positive. 



First order 

The isotropy at first order is described in the following Lemma. 

Lemma 10. The D2Q9 scheme is isotropic at first order if, and only if the 
properties (I1-I2) are satisfied. 

(11) At equilibrium, the energy is proportional to the density (the propor- 
tionality factor is relative to the sound velocity). 

(12) Both moments m^^^ and m^^^ vanish at equilibrium. 

Properties (I1-I2) involve the following structure for the matrix E (there is 
no constraint on the matrix S at first order): 



E 



/E§X^ 

E§,\^ 

E^PyX 





E^X^ 
EllX^ 






\ 

E'^IX^ 



Remark 11. This result gives a rigorously justification for usual physical 
assumptions. 



Proof. Considering eqs. eq211coslsinl and eq211cos0sin2 yields to -E^ 
Elpxy = 0. Then the system of four equations 



y>xx 



eq222cos2sinl = 

eq222coslsin0 = 

eq222cos0sinl = 

eq222cos0sin0 = 

gives that the four unknowns E^l^ , E^l^ , E^l^ 

To solve all of the remained equations, it is sufficient to consider the 
equations eq212cos0sin0 and eq322cos0sinl in the unknowns and 



0, 
0, 
0, 
0, 

and E^^^y vanish. 



Es^ : more precisely the equilibrium energy is characterized by E^ = El 
0. 



□ 



Second order 

For the D2Q9 scheme, the isotropy at second order can be characterized 
by the following Lemma: 

Lemma 12. The D2Q9 scheme is isotropic at second order if and only if 
the properties (I1-I4) are satisfied. 

(13) At equilibrium, the heat flux is proportional to the momentum, that is 
m^x = E^l^'^Qx and m^y = E^lX^qy. 

(14) The proportionality factor E^'l is related to the relaxation times by the 



relation E'^^ 



-(4cj, 



'^<Pxx)/{^'Pxx + 2cr, 
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Properties (I1-I4) involve the following structure for the matrix E (there is 
no constraint on the matrix S at second order): 



E = 



(ElX^ 





^ 


E'e,\' 


EllX^ 


EllX^ 





cA2 






















I 





) 



with c = -(4o-^,^ - o-^,J/(fT<^,., + 2cT^,J. 

Proof. Considering first eq. eq2221 and more precisely the system 

J eq2221cos0sin0 = 0, 
\ eqSlllcosOsinO = 0, 

we obtain that Ef^^, E^^ vanish. Then the system 

J eq2222cos0sin2 = 0, 
\ eq3112cos0sin2 = 0, 

imphes both results = E'^^ and i?^^ = —E^^.. Using eq. eq3113cos0siii2 
yields to E!^^ = and eq. eq2233cos2sin2 gives the characterization on 
-E^^. Finally, we prove that the previous relations are sufficient to impose 
the isotropy at second order. □ 

Remark 13. The relation E^l = -{^o^^y - '^^p^xMi.'^^^:^ + '^'^<fixy) ^■^ exactly 
the relation (4I) of [TT] except that the roles of the relaxation times a^^^ 
and (Tip^y have to be exchanged. 

Third order 

For the D2Q9 scheme, the isotropy at third order can be characterized by 
the following Lemma: 

Lemma 14. The D2Q9 scheme is isotropic at third order if and only if the 
whole of the properties (I1-I5) and one of the properties (Ig-Is) are satisfied. 

(I5) The relaxation times a^^^ and a^^^ (relating to second order moments) 
are the same, so that c = — 1. 

(Ig) At equilibrium, the square of the energy is proportional to the density: 
that is Eel = o,nd E^l = 0, and is linked with the energy through the 
relation 2E§2 +i + 3E§ = 0. 

(I7) At equilibrium, the square of the energy is proportional to the density: 
that is Eel — ^'''^d Ee^ = 0, and both relaxation times related to the 
heat flux satisfy a^^ = a^^ = 1/(12(7^^^). 

(Ig) Both relaxation times related to the heat flux satisfy a^^ = a^p^ = 
l/(12(7(p^^), and both viscosities are linked by ae = (^ip^^- 
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Properties (Ii"l5)+(l6) involve the following structure of the matrices E and 
S: 

( E^X'^ \ 

" ' I 

E= -A2 

-A2 



/ 



V 

5 = Diag {s. 



e ) ) Sip^ , Sipy , Sip^^ , S^p^^ 



Properties (Ii-l5)+(l7) involve the following structure of the matrices E and 



E 



S = Diag Se,Se2,3 



(E^X^ 





\ 


E'e.X^ 











-A2 











-A2 











V 





/ 








3 ~ Sipxx 




^^xx J 



Properties (Ii-l5)+(l8) involve the following structure of the matrices E and 



E 



(E^eX^ 
E§,X' 




V 



S = Diag s^^^,Se2,3: 






\ 


E!-,x^ 


i?gA3 


-A2 








-A2 











/ 


^fxx q2 


~ ^ifixx 


^fxx "J 


~ ^fxx 



X 1 °<fxx 



Proof. We first consider eq. eqllll2cos2sin2 and we establish property 
(le). Then, multiple choices appear considering eqs. eq32223, eq22221 and 
eq21112. In order to give more details about these choices we have to define 
seven relations: 



(la) 
(lb) 
(Ic) 



Etl 

0-e = 

TPQx 
^S2 



0, 



Vxx': 



0, 



(Id) 

(le) 
(If) 

(Ig) 



4 + 2E§2 + 3E§ = 0, 

-1 + 4o-<^..(2crv'. + (^Vy) 
0, 



0, 



'Vx 



cr, 



Vy- 



The investigation of these equations involves some of these seven relations. 
More precisely, we have: 

- eq. eq32223coslsin0 yields to a dichotomy between (la) and (Ic 

- eq. eq32223cos0sinl yields to a dichotomy between (lb) and (Ic 

- eq. eq32223cos5sin0 yields to a dichotomy between (la) and (Id) 

- eq. eq32223cos4siiil yields to a dichotomy between (Ic) and (Ij) 



Ic), 
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- eq. eq22221coslsinl yields to a dichotomy between (le) and (If), 

- eq. eq21112cos0sin2 yields to a dichotomy between (le) and (Ig). 

That gives a priori 2^ = 256 possibilities. However it is straightforward 
that only three cases remain: 

- properties (la), (Ic) and (le) are true, 

- properties (la), (Ic), (If) and (Ig) are true, 

- properties (lb), (Id) and (If) are true. 

Finally, each of these three cases solves all of the reminded equations and 
we get isotropy at the third order. □ 

Remark 15. 

- The properties (Is-Ig) are already used in the literature (the equality 2E§2 + 
4+3£'e = is a generalization of the classical choice E§ = —2 and = 1 
so it is really interesting to obtain a justification of this choice by the study 
of the isotropy. 

- The equality s^^ = 3(2 — s^^^)/{3 — s<^^^) is also proposed in [H]. 

- To our knowledge, the sets of parameters given by properties (Ii-l5)-(l7) 
and (Ii-l5)-(l8) are new though the second one is less interesting for acoustic 
applications because of the link between both viscosities. 

Fourth order 

For the D2Q9 scheme, the isotropy at fourth order can be characterized 
by the following Lemma: 

Lemma 16. The D2Q9 scheme is isotropic at fourth order if and only if 
the whole of the properties (If-Ig), (Ig-In) and one of the two properties 
(Ii2^Ii3) are satisfied. 

(I9) The relaxations times related to odd moments satisfy a^^ = = 

(110) At equilibrium, the square of the kinetic energy is proportional to the 
density E2 = E^^p}^ . 

(111) At equilibrium, the kinetic energy and its square are linked by 2E§2 + 
4. + 3E§ = 0. 

(112) The relaxations times related to even moments satisfy as = 0-^2 = (^ip^x- 

(113) The sound velocity is imposed through the equality 2 + 3E§ = 0. 
These conditions involve the following structure for the matrices E and S: 



E 



2 








5 = Diag s„s,,,6- ^, 

6 - s^^^ 






\ 



























/ 













where either 3^2 = Se or E^ = —2/3. 
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Proof. Since there are three sets of coefficients giving third order isotropy, 
the proof is divided into three cases. 

The cases with (I1-I5) and either (I7) or (Ig) are forbidden because they 
both imply the condition o'^^^ = 0, using eq. eqllllllcos2sin2. 

So we have to assume (Ii-Ie) in order to impose isotropy at third order. 
These assumptions imply E'lj = EeJ^ = using eqs. eqlll222cos5sin0 
and eqlll222cos4sinl and a^y = a^^ by eq. eq212222cos5sinl. Then 
eq. eq311113cos2sin2 involves property (I9) and eq. eqll2221coslsinl 
involves property (In). Finally the dichotomy between (I12) and (I13) comes 
from eq. eq211112cos2sin2. □ 

Remark 17. 

- This Lemma is the last step to prove that the heat flux and the momentum 
are collinear: the property of isotropy at fourth order imposes the equality 
between both relaxations times s^p^ and s^y. Moreover, it links these two 
relaxation times with the other one s^^^ . 

- The property (In) is involved in a particular case in [11] (E§ = —2, 
E§2 = but without the constraint on the relaxation times given in (I9). 

- The case (Ii-l5)-(l9-Iii)-(Ii3) insures isotropy at fourth order. However, 
the combination of relation E§ = 6cg— 4 and property (I13) yields to Cq = 5/9, 
and this value for the sound velocity involves that the bulk viscosity disap- 
pears by using eq. (fTTjl [6]. 

- The isotropy at fourth order is quite restrictive: both relaxation times linked 
to the viscosities (Sg and s^p^^) must be equal. 

- Note that these results generalize those obtained in ^ for quartic param- 
eters. More precisely, the values a^^ = -v/3/3 and a^^^ = \/3/6 proposed in 
[9] are compatible with property (I9). 

Fifth order 

For the D2Q9 scheme, we have the following Lemma: 

Lemma 18. The D2Q9 scheme is never isotropic at fifth order. 

Proof. In both cases (Ii-l5)-(l9-Ii2) and (Ii-l5)-(l9-Iii)-(Ii3), Proposition 2] 
gives equations that can not be solve for every rotation of the frame: for 
example eqllllll2 is of type ccos^sin^, where c is a given real constant 
independent of the parameters (namely, c do not vanish) . □ 

3.3. D2Q13 scheme 

In this Subsection, we give the main result on the D2Q13 scheme in Propo- 
sition [T21 The proof of this Proposition is then subdivided into two Lemmas 
in order to detail the methodology due to the definition of the isotropy. 

Proposition 19. Let Li{r) be the lack of isotropy for the D2Q13 scheme at 
fourth order for the rotation r, then we get the following propositions. 

• The scheme is isotropic at first order, that is Li{r) = C'(At^),Vr E 
SO^m, iffEr, El\ E^^^, E$l^, EX., EP^^y, E^^, and E^^ are 
zero. 

• The scheme is isotropic at second order, that is Li{r) = 0{/S.t^),\/ r E 
S'02(K), iff the four following properties are satisfied: 

— it is isotropic at first order. 



16 ORSAY, ORSAY.CNAM, AND ORSAY 

- E!J,^ , ElJ,y , E^e2 ) "'iT'd E^e2 a'^e zero, 

- E'j,l and E^l are equal to (a^^^ - 4ct^,J/(3o-^^^ + cr<^^J, 

_ rfs -ply T^Hx „„J T^'iy 

-^xe2 — ''--"y£2! -"^fx — fy' ''-'xe2 — -^J/£2 • 

The proof of Proposition [19] is detailed below thanks to two Lemmas: 
one for each order. The methodology is the same as the one for the D2Q9 
scheme. 

First order 

For the D2Q13 scheme, the isotropy at first order is characterized by the 
following Lemma: 

Lemma 20. The D2Q13 scheme is isotropic at the first order if and only if 
the properties (I'^-l2) are satisfied. 

iVi) At equilibrium, the energy is proportional to the density (the propor- 
tionality factor is relative to the sound velocity). 

(I2) At equilibrium, both moments m^^^ and m^^^ are zero. 

Both properties (I^-l2) implies the following structure of the matrix E (there 
is no constraint on the matrix S at first order): 



( E§X^ 





\ 






















Tpqx \2 


EllX^ 


E^y^^ 






Exe2X 




jply \A 
^XS2^ 


Eys2X 


Tpgx x 4 


TpQy \4 
rjye2'\ 


E§2X^ 


Er2X^ 


EtlX^ 




EllX^ 


EtlX^ 


\ExxeX 


ExxeX? 


Et^xeX? ) 



Remark 21. The first order isotropy conditions are the same for the D2Q9 
and the D2Q13 schemes. 

Second order 

For the D2Q13 scheme, the isotropy at second order is described in the 
following Lemma: 

Lemma 22. The D2Q13 scheme is isotropic at second order if and only if 
the properties (I'l-Is) are satisfied. 

(13) At equilibrium the heat flux is a rotation- dilatation of the momentum, 
more precisely m^^ = E^lX'^q^ + E^^X'^qy and m^^ = -E'^^X^q^ + 

Et>^hy. 

(14) At equilibrium the moment of order five (mxeji "^?/e2) ^-^ ^ rotation- 
dilatation of the momentum, more precisely m^^g = Exe2X'^qx+Ex%2X'^qy 
and mys2 = -El^X^q^ + E^e2X'^qy 
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(I5) The equilibrium states are related to the relaxation times by the rela- 
tions Exe2 — 0- and Exe2 — with 

(14) 

Properties (l'^-l5) implies the following structure of the matrix E (there is 
no constraint the matrix S at second order): 



E = 



( E§X^ 





\ 























EtX' 


Et^^ 





-EtX' 


E%X' 





aX* 


bX' 





-bX^ 


aA^ 


E§2X' 


E^X^ 


EllX^ 




EllX^ 


EllX^ 


\ExxeX 


ExxeX? 


ExxeX? J 



where a and b are given in (|14p . 

Remark 23. The properties of isotropy at second order for the D2Q9 and 
D2Q13 schemes are of the same type, though they are less constraining for 
the D2Q13 scheme (the coefficient E^^^ does not have to be zero for instance). 

Third order 

Because of the very high number of cases offered to guarantee the isotropy 
at third order (to our knowledge, not less than 17 different cases have to be 
investigated), only some sufficient conditions are given in this paper. First 
of all, we already know that the isotropy at third order implies properties 

(I7) At equilibrium, the heat flux is proportional to the momentum: that is 
m^x = E%X^qx and m^^ = Ef^X^qy. 

(Ig) At equilibrium the moment of order five (m^^^j , m^^j ) is proportional to 
the momentum: that is m^ea = Exe2X^qx and m^y = Ext2X'^qy. 

(Ig) The proportional coefficients and Ex%2 CLf^ linked by the relation 
Ell = -(21/8)^^: - 65/24. 

Assuming properties (I'^-Ig), we give some example of sets of coefficients 
that make the D2Q13 scheme isotropic in the Annex. 

4. Numerical results for the D2Q9 scheme 

In this Section, we present some preliminary numerical results in order 
to appreciate the lack of isotropy order after order. Since it is not easy to 
observe the lack of isotropy in the sense of Definition[5l we have to investigate 
an other way to make explicit the isotropy error. 

Then, it seems natural to represent the evolution of a radial function after 
a few time steps. Let r 1— )■ p{r, 6) be the density solution of the linear lattice 
Boltzmann D2Q9 scheme described in Section 13.21 and initialized with the 
Gaussian for the first moment p{x, y) = exp(— lOx^ — lOy^) and with zero for 
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the others. Indeed the most the scheme is isotropic, the most the behavior 
of r I— )• p{r, 6) is the same for all of the angle 9. We consider a D2Q9 scheme 
with 100 X 100 space meshing (in our case Ax = 0.02 and A = 1), our 
purpose is to make explicit the isotropy error at t = 12At. Thus, we plot 
the density r 1— )■ p{r, 6) where 9 is fixed: 

• pq := r ^ p{r, 9 = 0) (drawn with V), 

• /Oi := r I— )• p{r,9 = ^) (drawn with A), 

• p| := r I— 7- p{r,9 = j) (drawn with o), 

• Parctan{i/2) ■= r ^ p{r,9 = arctan(l/2)) (drawn with <). 

For each density, we only have the values on the nodes of the mesh: namely 
for each angle 9, the vectors of the abscissa of {r^pe) are not equal. Since 
we have to compare these densities, we interpolate the value with a spline 
method of order 5 so that the interpolation error does not interfere. 

We first consider the following parameters for which the scheme is isotropic 
at first order but not at second order: 



E 



/ 


-2 





o\ 




6 













-2 













-2 













V 








0/ 



0.5s, 



and S = diae 



where C = 1.84 10-^ p = 2/3C, s, = 3C/(AAx) = 1.9977944349438221, 
s^^ = 1.3, s^^^ = 3^/(AAx) = 1.9985290825952098. 

The different curves pg are plotted in Figured) The isotropy of the scheme 

is weak in the sense that, for instance, the plot of po — pi in Figure [2] gives 

4 

an error of size 6.5E-4. We remark that the curves of po and piL are the 
same: the scheme is invariant by a rotation of 7r/2. 




Figure 1. Isotropy 
at first order 




Figure 2. Isotropy 
error at first order 



We then consider the following parameters, with no change on the relax- 
ations times, for which the scheme is isotropic at second order but not at 
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third order: 



E 



/ 


-2 





o\ 




6 













-1 













-1 













V 








0/ 



and S = diag 



'''fix 
'^^xx 



\ ^^xx / 



In Figure El we observe that this choice of parameters improves the 
isotropy in the sense that the different curves are more similar but the 
difference between po ^-nd pn is still of order 4.5E-4 on Figure HI 





Figure 3. Isotropy 
at second order 



Figure 4. Isotropy 
error at second order 



We next consider the following parameters, with no change on the re- 
laxations times, for which the scheme is isotropic at third order but not at 
fourth order: 



E 



( 


-2 





o\ 




1 













-1 













-1 
























and S = diag 



0.5s, 

^fx 
^fx 
^^xx 
\ ^ipxx / 



The results are given in Figure [5] is better than the one given for isotropy 
only at the second order: the plot of po ~ P- in Figure [6] gives an error of 
size 4.2E-6. 





Figure 5. Isotropy 
at third order 



Figure 6. Isotropy 
error at third order 
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And finally, we consider the following parameters for which the scheme is 
isotropic at fourth order: 



E 



(-2 





o\ 


1 











-1 











-1 











Vo 





0^ 


-5 _ 




= 3 



and S = diag 



(Se\ 



Se 
\SeJ 



3C/(AAx) = 1.9977944349438221, s, 



iPx 



where C = 1.8410 

6.0(2.0 - s^,J/(6.0 - s^^J = 0.0022055650561781941. 

The result is given in Figure [7] and the plot of po — in Figure [8] gives 
an error of size 5.5E-7. 





Figure 7. Isotropy 
at fourth order 



Figure 8. Isotropy 
at fourth order 



CONGLUSION 

We used a general methodology that guarantees isotropy of a lattice Boltz- 
mann scheme for a given order. This approach is based on the expansion of 
the equivalent PDEs at an arbitrary order and the invariance by the group 
operation We have detailed all the possible cases for the basic scheme 
D2Q9 applied to a linearized fluid mechanics. Results (up to the second or- 
der) have also been proposed for the D2Q13 scheme. Using "inappropriate" 
choices of parameters in the D2Q9 scheme, elementary test cases hightlight 
the lack of isotropy at various orders. This work can be extended without 
conceptional difficulties for 3D lattice Boltzmann schemes and it will be 
done in a forthcoming contribution. 



Annex 

In this annex, we propose some details concerning the isotropy at third or- 
der for the D2Q13 scheme. This work is in progress: we express the different 
cases that have to be study. 

In fact, knowing properties (I'^^-Ig), eqs. eqllllS and eq22233 give mul- 
tiple choices between the thirteen following relations: 
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(la) 


Et = -3, 


(If2) 
(Ig) 


(ID 




(1^1 ) 


First value for E^^, 


(Ifs) 


(Id) 


aye2 = 1/12<T<^^^, 


(i;.2) 


(i;) 




(I^s) 


(Ifl) 


First value for Ej^ , 


(4) 


(Ic2) 


Second value for Ee^, 





Second value for E^^, 

^{(^Vxx + '^f>xy){'^xe2 + ^ye2) = 1' 
Third value for E^^, 

^^fxy {^xe2 ~^'^ys2 

) + ^(^<fixxi(^ye2+^cr 
Third value for E^^, 
Fourth value for Es^. 



The investigation of these equations involves some of these thirteen rela- 
tions. More precisely, we have: 

- eq. eqllllScoslsinl yields to a dichotomy between (l'^) and (![,), 

- eq. eq32223cos0sin0 yields to a dichotomy between (I^^^) and 

- eq. eq32223cos0sinl yields to a dichotomy between (Ig) and (If^), 

- eq. eq32223coslsin0 yields to a dichotomy between (l'f.2) and (Ig), 

- eq. eq32223cos4sinl yields to a dichotomy between (Ijg) and (Igi), 

- eq. eq32223cos2siiil yields to a dichotomy between (Ijg) and (Ig2), 

- eq. eq32223cos3sin0 yields to a dichotomy between (1^3) and (Ig2); 

- eq. eq32223cos5sin0 yields to a dichotomy between (Igi) and (1^4). 

That gives a priori 2^ = 256 possibilities. However preliminary calcula- 
tions yield that only seventeen cases remain: in order to have isotropy of 
third order, it is necessary to satisfy at least one of these sets of properties. 

1) (I'J, (4), (4), (Id), (le), (4) and (4) 

2) {!',), (I^), (I^), (i;.,) and (4) 

3) (I'a), (4), (4), (4), (ID, (Ifl), (4) and (4) 

4) (I'J, (4), (4), (ID, (4), (4) and (4) 

5) (ID, (4), (4), (ID, (4), (4) and (4) 

6) (I'J, (i:,J, (4), (4), (IJ, (4) and (4) 

7) (I'J, (i;,J, (4), (IJ, (4) and (4) 

8) (I'J, (4), (4), (i:,), (4) and (4) 

9) (IJ, (4), (4), (4), (4), (4), (4) and (4) 

10) (IJ, (4), (4), (4), (4), (4) and (4) 

11) (IJ, (4), (4), (4), (4), (4) and (4) 

12) (IJ, (4), (4), (4), (4) and (4) 

13) (IJ, (i:,J, (4), (4), (4), (4), (4) and (4) 

14) (IJ, (4), (4), (4), (4), (4) and (4) 

15) (IJ, (4), (4), (4), (4), (4) and (4) 

16) (I'J, (4), (4), (4), (IJ, (4), (4) and (4) 

17) (IJ, (IJ), (4), (4), (IJ, (4) and (4) 

The study of these cases is in progress and we cannot then give a full 
characterization of it. However we know some examples of matrices E and 
S that involve isotropy at third order. We propose here two cases given by 
properties (Ij, (4), (4), (Ij, (Ij, (4) and (4). 
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E 



274 67 
^ ~ 462 








° ^ 























-3A2 











-3A2 





31 4 

¥^ 










-A^ 

EtlX^ 




EllX^ 


-EP, Ei] 
' 3003 J 


462 ' 


"462^^^^ 






ExxeX? ) 



diag (Sip^^ , Sip^^ , S^^j^ , S^p^ , S^p^ , Sip^, S^p^ , S^^ , S^g , Sxxe) ) 



where s^, = 3- — 



or for instance 



E 



ElX^ 









3234 361 



E'e 



X' 



13 26 
/1681 307 ) 

V 39 156 




S = diag(^S£,s^^^,s<^^^,2 




















-3A2 








-3A2 


31 4 






6 

















/ 



where 5 



-As 



2-s. 



) Sipy 1 Sx£2 1 ^xe2 1 ^£2 ) '^£3 ' ^XX£ ^ 



xe2 



Remark 24. Since the D2Q13 scheme takes into account the velocities of 
the D2Q9 and 4 additional ones, it can be read as a generalization of the 
D2Q9 scheme. In order to illustrate this remark, we propose here a set 
of coefficients that gives isotropy at third order for both D2Q9 and D2Q13 
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schemes: 



E 



E§X^ 


















































-A2 







12 




-i-A' 

12 






3234 361 
13 26 % 








ri681 307 A 
I 39 156 / 


A6 









V 











/ 



S — diag{s^p^.^, s^p^^, ^^pxxi '•Px-, Sip.^., s^p^, s^p^, s^^i s^^, oxxe 

^2 s<n^^ 



where s, 
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